Invariant submanifold for series arrays of Josephson junctions 
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We study the nonlinear dynamics of series arrays of Josephson junctions in the large- N limit, where 
N is the number of junctions in the array. The junctions are assumed to be identical, overdamped, 
driven by a constant bias current and globally coupled through a common load. Previous simulations 
of such arrays revealed that their dynamics are remarkably simple, hinting at the presence of some 
hidden symmetry or other structure. These observations were later explained by the discovery of 
N — 3 constants of motion, each choice of which confines the resulting flow in phase space to a 
low-dimensional invariant manifold. Here we show that the dimensionality can be reduced further 
by restricting attention to a special family of states recently identified by Ott and Antonsen. In 
geometric terms, the Ott- Antonsen ansatz corresponds to an invariant submanifold of dimension one 
less than that found earlier. We derive and analyze the flow on this submanifold for two special cases: 
an array with purely resistive loading and another with resistive-inductive-capacitive loading. Our 
results recover (and in some instances improve) earlier findings based on linearization arguments. 



Josephson junctions are superconducting de- 
vices with many practical applications, ranging 
from voltage standards to ultrasensitive detec- 
tors. From a mathematical perspective, their 
nonlinear dynamics are fascinating, especially 
when many junctions are coupled together in an 
array. For about the past twenty years, theorists 
have been intrigued by the strange collective be- 
havior seen in numerical experiments on arrays of 
identical junctions in series. The behavior began 
to make sense when it was eventually realized that 
despite the presence of dissipation in the under- 
lying circuits, the equations possess an enormous 
number of constants of motion. These constants 
restrict the dynamics to low-dimensional mani- 
folds in phase space. In this paper, we show that 
in certain cases the reduced dimensionality can 
be even more severe than previously realized. By 
making use of an ansatz recently introduced by 
Ott and Antonsen, we demonstrate for example 
that a resistively loaded array can behave exactly 
as if its phase space were two-dimensional, even 
when the array consists of infinitely many junc- 
tions. 



I. INTRODUCTION 

Forty years ago, Winfree pioneered the study of syn- 
chronization in large populations of coupled limit-cycle 
oscillators [1]. Since then, the field has expanded con- 
siderably, thanks in large part to Kuramoto's elegant 
reformulation [2, 3] of Winfree's intuitive model. Both 
of their models were originally motivated by biologi- 
cal phenomena [4] such as the alpha rhythm of brain 
waves [5, 6, 7], the collective firing of cardiac pacemaker 
cells [8, 9, 10], the coordinated flashing of southeast Asian 
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fireflies [11, 12, 13, 14, 15], menstrual synchrony among 
close female friends [16, 17, 18], and glycolytic oscillations 
in yeast populations [19, 20, 21, 22, 23, 24]. However, 
the techniques developed in analyzing these systems soon 
proved relevant to physical problems, including the dy- 
namics of laser arrays [25, 26, 27, 28, 29, 30, 31], charge- 
density waves [32, 33, 34] and coherence among sites of 
electrochemical dissolution [35, 36, 37]. 

For simplicity, all the individual oscillators in such 
models have often been assumed to be coupled equally 
strongly to all the others, a form of interaction known 
variously as global, infinite-range, or mean-field coupling. 
While this form of coupling is a crude approximation in 
most cases, series arrays of Josephson junctions consti- 
tute a notable exception. In these systems, exact global 
coupling between the Josephson junctions emerges nat- 
urally from Kirchhoff's laws and the physical properties 
of weakly coupled superconductors [38, 39]. 

In the early 1990s, numerical simulations of Joseph- 
son junction arrays revealed that they were prone to a 
large degree of neutral stability [38, 39, 40, 41, 42]. This 
peculiar phenomenon was first seen in a simple array in 
which N identical junctions were connected in series and 
coupled by a resistor in parallel with all of them. The 
numerics suggested that the system's trajectories were 
always trapped on two-dimensional tori, no matter how 
many junctions were included in the array [38]. Later 
studies showed that arrays of identical junctions coupled 
through other kinds of loads displayed a similarly non- 
generic form of behavior: all but four of the Floquet mul- 
tipliers for a certain periodic state known as a splay state 
appeared to lie on the complex unit circle and remain 
there as the system parameters were varied [39, 41]. 

These and other puzzling observations were partially 
explained by the subsequent discovery that the equa- 
tions of motion could be reduced to a dynamical sys- 
tem of much lower dimension [43, 44, 45]. Specifically, 
Watanabe and Strogatz [43, 44] found a time-dependent 
trigonometric transformation that expressed the junction 
phases in terms of N constant phases and three collective 
time-dependent variables, each obeying a suitable differ- 
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cntial equation. The transformation took the form 
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for Josephson junctions j = 1, . . . , N. Here, the variables 
(j)j{t) denote the junction phases, while the constants 0j 
denote the fixed phases on which the transformation op- 
erates, and 7(f) and <d(t) denote the collective vari- 
ables. 

In the limit of infinitely many Josephson junctions, an- 
other result could be extracted from (1). The transfor- 
mation maps old phases 6j to new phases <f>j in such 
a way that a uniform distribution in 8, with phases 
spread evenly around the circle, is transformed into a 
non-uniform distribution in cf>, with phases symmetrically 
clumped about some mean phase and distributed accord- 
ing to a Poisson kernel. This result implies that the space 
of Poisson kernels is dynamically invariant. The reason- 
ing is as follows: if the phase distribution takes the form 
of a Poisson kernel at any time, it can be viewed as having 
arisen (via the transformation) from an initially uniform 
phase distribution and hence will remain distributed as 
a Poisson kernel for all time. 

Recently, Ott and Antonsen [46] showed by explicit 
calculation that the submanifold of Poisson- kernel phase 
distributions — henceforth called the Poisson submani- 
fold — is invariant for a much wider class of models, 
including the Kuramoto model and other systems in 
which the oscillators are non-identical. For the Ku- 
ramoto model in which the oscillators have Lorentzian- 
distributed natural frequencies, Ott and Antonsen [46] 
derived an exact differential equation for the evolution of 
the complex order parameter (the centroid of the phase 
distribution around the unit circle). As it happens, the 
amplitude and phase of the order parameter completely 
characterize the Poisson kernel. Hence, by extracting 
the dynamics of the order parameter, Ott and Anton- 
sen [46] simultaneously unveiled the dynamics on the 
Poisson submanifold. 

Our goal in this paper is to use Ott and Antonsen's 
ansatz [46] to investigate the dynamics of series arrays 
of identical Josephson junctions. Before turning to this 
specific task, however, we first consider the scope of the 
ansatz itself. What algebraic form do the governing equa- 
tions need to have in order for the ansatz to work? In 
Sec. II, we pinpoint the family of equations that can be 
simplified in this way, a family that includes Josephson 
junction series arrays as a special case. 

Sections III and IV then apply the ansatz to two par- 
ticular arrays, one with a resistive load and another with 
an RLC load. The analysis illuminates the order param- 
eter dynamics on the entire Poisson submanifold. In this 
way, we clarify for the first time how the associated ar- 
rays behave, not just near their equilibria and periodic 
orbits, but also far from those special states. 

Finally, Sec. V discusses what the ansatz does — and 
does not — imply about the dynamics of the original 



Josephson junction arrays. Although the ansatz pro- 
vides powerful insights, it does not tell the whole story; 
it misses certain important dynamical states that lie off 
the submanifold of Poisson kernels. These more general 
states can be handled by considering (1) within the for- 
malism of Mobius transformations, as has recently been 
shown elsewhere [47, 48] . 



II. REDUCIBLE SYSTEMS 

The most extensively studied systems of phase oscilla- 
tors, from the Kuramoto model to Josephson junction ar- 
rays, involve purely sinusoidal interactions. This single- 
harmonic structure is the key to the success of the Ott- 
Antonscn ansatz [46], as the following calculations show. 

Consider a system of N identical phase oscillators gov- 
erned by 



4> 3 = fe-^ +g + fe- 
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for j = 1, . . . , N . Here / is any smooth, complex- valued, 
27r-periodic function of the phases (j>\, . . . , 4>n- The func- 
tion / is allowed to depend on time and any other auxil- 
iary state variables in the system (for example, the charge 
on a load capacitor or the current through a load resistor 
for the Josephson junction arrays discussed in Sec. IV). 
What is crucial, however, is that / must not depend on 
the oscillator index j; it must be the same function for 
all j. Likewise, the function g must be independent of j. 
Note that g has to be real-valued since <fij is real. 

Intuitively, the functions / and g can be regarded as 
common fields felt by all the oscillators. These fields 
might involve averages over all the phases, as in mod- 
els with mean-field coupling, but this is not necessary. 
In fact, (2) need not even have permutation symmetry. 
That is, the equations need not stay the same under an 
arbitrary interchange of indices, because the functions 
/ and g need not respect such a symmetry. The only 
requirement is that / and g must be the same for all j. 

The class of systems (2) was studied previously by 
Watanabe and Strogatz [44] and by Goebel [45], who 
showed that all equations of this form are solved by the 
transformation (1), where the evolution of 7 and O is 
governed by the forms of / and g. We now show that the 
newly discovered Ott- Antonsen ansatz also works on this 
same class of systems, suggesting some intimate relation- 
ship between the two reduction methods, as addressed 
by [47, 48]. 

The Ott- Antonsen ansatz is restricted to the infinite- A'' 
limit of (2). In this limit, one describes the system not in 
terms of the motion of individual oscillators, but rather 
in terms of the evolution of the phase density /?(</>, i), 
defined such that p((j),t)d(j) gives the fraction of phases 
that lie between <p and (f> + &<j) at time t. Then p satisfies 
the continuity equation 
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where the velocity field is 

v(<l>,t) = fe- i ++g + fe i +, (4) 

from (2). Here, in the case of infinite N, our assumptions 
about the coefficient functions / and g take the form that 
/ and g may depend on t but not </>. The time-dependence 
of / and g can arise either explicitly (through external 
forcing, say) or implicitly (through the time-dependence 
of the harmonics of p or any auxiliary state variables in 
the system). 

Next, following Ott and Antonsen [46], suppose p is of 
the form 

K<M) = ^{l + £>(i)V"* + a(t) n e-^)\ (5) 

^ n— 1 ' 

for some unknown function a that is independent of <j). 
Note that (5) is just an algebraic rearrangement of the 
usual form for the Poisson kernel: 



1 2?r 1 - 2r cos(0 - ifj) + r 2 w 

where r and tp arc defined via 

a = re- tlp . (7) 

In geometrical terms, the ansatz (5) defines a subman- 
ifold in the infinite-dimensional space of density func- 
tions p. This submanifold is two-dimensional and is 
parametrized by the complex number a (or equivalcntly, 
by the polar coordinates r and ip). 

We now show that the submanifold of Poisson kernels 
is invariant and calculate the flow on it. To do so, we 
substitute the velocity field (4) and the ansatz (5) into 
the continuity equation (3). After reindexing where ap- 
propriate, we obtain 

oo 

[d + i{f + ga + fa 2 )] £ na n - 1 e in<f ' + c.c. = (8) 

n=l 

where c.c. denotes the complex conjugate of the first term 
in (8). Note the seemingly miraculous coincidence here: 
the expression in brackets is a common factor for each 
term in the sum. Hence, although (8) constitutes an 
infinite set of amplitude equations, with one for each 
harmonic, all of them are satisfied simultaneously if the 
bracketed expression vanishes. This condition is also nec- 
essary, since 

°° i<j> 

a-^n(ae^= ^0. (9) 

71=1 ^ ' 

Thus (8) is satisfied for all <f> if and only if 

a + i{f + ga + fa 2 ) = 0. (10) 

It proves convenient to reexpress this result in terms of 
the complex order parameter z, defined as usual by the 
ccntroid of the phase distribution: 

Z (t)= r > P (<M)c¥. (ii) 

Jo 




FIG. 1: A series of Josephson junctions in parallel with a 
resistive load. /(, is the constant bias current, R is the load 
resistance, and Rj is the internal resistance of a single Joseph- 
son junction. 

Then by substituting (5) into (11) we find that z — a = 
re 1 ^ . Hence, z satisfies the Riccati equation 

z = i(f + gz + fz 2 ). (12) 

This equation gives the flow on the Poisson submanifold. 
When / and g are functions of z and t alone, as in the 
case that they can be expressed in terms of the harmonics 
of p, (12) constitutes a closed two-dimensional system. 

III. RESISTIVELY-LOADED CIRCUIT 

We now apply (12) to obtain the reduced dynamics 
of the circuit shown in Fig. 1. This circuit consists of A" 
Josephson junctions wired in series and placed in parallel 
with a resistive load R and a constant current source lb- 
Let I c denote the critical current of each junction and 
let Rj denote its internal resistance. For simplicity, the 
junctions are assumed to be heavily overdamped so that 
we can neglect their internal capacitance. 

As demonstrated by Tsang et al. [38], the time- 
dependent dynamics of this circuit can be converted via 
Kirchhoff's laws and the physical properties of Josephson 
junctions to the dimcnsionless system: 

1 N 

<pj = + a cos <f>j + — ^ cos <f>k (13) 
fc=i 

for j = 1, . . . , N. Here, <pj = 6j — ir/2, where 8j is the 
phase difference across the jth oscillating Josephson junc- 
tion. We write the system in these unusual variables to 
highlight its reversibility symmetry; the equations stay 
the same if we change <f>j — > —<f>j and t — > —t. The di- 
mcnsionless groups a and CI in (13) are given in terms 
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of the original circuit parameters by a = —RJ (NRj) — 1 
and Q = I b R/(NI c Rj). Taking the limit N -> oo of (13), 
we obtain the velocity field 



If 6 7^ 0, ±1, (17) has four fixed points: 



= ~e + (Q + rcos-0) + -e 



(14) 



where r and r/> denote the amplitude and phase of the 
complex order parameter z, as before. Equation (14) has 
the form of (4), so by (12), the order parameter dynamics 



— + [il + rcosipjz + — z 



(15) 



Taking the real and imaginary parts of (15) yields the 
two-dimensional nonlinear system: 



■ l + b ( 2 U ■ I 

r = — - — (r — 1) smip 

tp = Q i— (r + r -1 ) cos ip + r cos ip 



(16) 



where we have introduced the parameter 6 = —a — 1 = 
R/(NRj), which simplifies the parameter space and fixed 
points of (16). Note that b > for any real circuit, al- 
though in what follows we will also allow negative values 
of b, since our main interest is with (13) as a dynamical 
system, not as a model of a real device. 



A. Fixed points 



x* = Q/b; y* = ^l-fl 2 /b 2 
x* = n/6; y* = -^/l-fl 2 /b 2 



. -n + yn 2 -b 2 + i , 

x = ; y = 

1 - b ' y 



_ n _ ^2 _ b 2 + 1 
x = — ;y =0 



(18a) 
(18b) 

(18c) 
(18d) 



The points (18a) and (18b) lie on the boundary r = 1 
of the unit disk and represent synchronized rest states: 
equilibrium states in which all the Josephson junctions 
are perfectly in phase and do not oscillate. In this case, 
all the individual phases <pj equal the same constant, and 
hence equal the constant phase ip of the centroid. Physi- 
cally, such states would be superconducting. The source 
current tunnels through each of the N junctions without 
developing any voltage across the load. 

In contrast, the fixed points (18c) and (18d) have r < 1 
and lie inside the unit disk, meaning that the Josephson 
junctions are not all in phase. This type of fixed point is 
known as a splay state [38, 39, 40, 41, 42]. It represents 
a periodic collective state in which the junctions oscillate 
out of phase, but in such a highly organized way that the 
overall phase distribution remains stationary. In particu- 
lar, the macroscopic order parameters r and ip stay con- 
stant even though individual junctions change their state 
non-uniformly, hesitating at some phases and accelerat- 
ing at others. The stationary distribution of phases takes 
the form of a Poisson kernel, as expected. 



This section analyzes the fixed points of (16) with re- 
spect to their dependence on the parameters b and O. 
There are many cases so we organize the calculations as 
a series of simple claims. Readers who prefer to bypass 
the algebra can skip ahead to Sec. HID which summa- 
rizes the findings. 

First we convert (16) to Cartesian coordinates. This 
removes the coordinate singularity at r = and permits 
simpler proofs of where the fixed points of (16) exist in 
the b-fl parameter space. Let x = r cos ip and y = r sin ip. 
Then (16) becomes 

x = bxy — fly 

. 1-6 2 _ 1 + 6, 1 + 6 (17) 

y = x + ilx H y" 

y 2 2 2 

Note that (17) remains identical under the transforma- 
tion Q — > — Q, x — > —x, which represents the symmetry 
of the resistively-loaded circuit in Fig. 1 under a sign re- 
versal of lb and reflection of the coordinate system on 
which phase synchrony is measured. Hence, the number 
and stability of the fixed points remains unchanged for 
each point in 6-f2 space reflected across the 6-axis. With- 
out loss of generality, we therefore consider only positive 
values of from here on. 



B. Partitioning the parameter space 

We now determine where in the 6-f2 parameter space 
the fixed points exist. By inspection, the synchronized 
rest states (18a) and (18b) exist if and only if 17 < |6|. 
Meanwhile, for the splay states (18c) and (18d) to exist, 
x* must be real (i.e. ft 2 - b 2 + 1 > 0) and < 1. The 
condition |x*| < 1 places additional restrictions on the 6 
and fl values at which (18c) and (18d) exist. We now 
derive these restrictions explicitly, with the key results 
summarized in Table I. 

Claim: On 0, 6 > 0, (18c) exists if and only if fl > b. 
Proof: Let h(b, ft) = tt 2 - b 2 + 1. Then by algebraic 
rearrangement , 

^ 1 ^ n ^^ n2 - b2 + b - ( 19 ) 

When the discriminant of (18c) and (18d) is nonnega- 
tive, fl 2 -b 2 + b > 6-1. Furthermore, 6-1 > for 6 > 1, 
and 6(1 — 6) > for 6 on [0, 1], so both sides of the second 
inequality in (19) arc nonnegative for all 6 > 0. Thus, 
we can square both sides of this inequality to remove the 
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TABLE I: Existence of the fixed points for S7 > 0. 



fixed point 


conditions 


(18a),(18b) 


n < \b\ 


(18c) 


Q, 2 - b 2 + 1 > and > b 


(18d) 


Q, 2 - b 2 + 1 > and fi < -6 



TABLE II: Linearization of (17) at (18a) for Q > 0. 



interval of b signs of A, t+ 


stability 
classification 


(-oo,min{-l,-0}) A > 0; r+ < 
(-l,min{-i/ 2 . -0}) A<0;r+<0 
(-V2,-fi) A<0;r+>0 
{Q,+oo) A>0;r+>0 


stable node 
saddle point 
saddle point 
unstable node 


TABLE III: Linearization of (17) at (18c), (18d) for fl > 0. 


fixed point interval of b sign of A± 


stability 
classification 


(18c) (-vo 2 + 1, fi) A+>0 


center 


(18d) (-Vtt 2 + 1, A_ < 
min{ — 1, — Q,}) 
(-1,-0) A_ > 


saddle point 
center 



square root: 

n 2 h > {n 2 - b 2 + b) 2 n > |6|. (20) 

We can also see that (20) implies (19) when ft > 0, 
because there the larger quantity of the first inequality 
in (20) is nonnegative and remains so upon removing 
the square in going from (20) to (19). Hence, the claim 
is proved. □ 



2 



x 



-1 





//r 2 -4A = 


b<-Q\ 





-1 -0.5 0.5 1 

A 

FIG. 2: A and r of the linearization of (18a) plotted para- 
metrically as a function of b for Q = 

Proof: Again by algebra, 




< 1 & nVh < -n 2 + b 2 -b. (22) 



We can square both sides of (22) to obtain: 

n 2 h < {fl 2 - b 2 +b) 2 n < \b\. (23) 

In the reverse direction, removing the square from the 
first inequality of (23) requires that f2 2 — b 2 + b < 0. 
However, 6(1 - b) > for b on [0, 1] and ft 2 - b 2 + b < 
implies h < for b > 1, so this inequality is not satisfied 
for b > 0. Nevertheless, O < -6 implies ft 2 - b 2 + b < 
for b < 0, so (18d) exists if and only if ft, > and il < —b. 
□ 



Claim: On il, —b > 0, (18c) exists if and only if ft, > 0. 
Proof: In the forward direction, we again obtain (19), 
which we split into two cases: 

flVh > \Q 2 - b 2 + b\ ^ fl > -b, (21a) 

nVh < \n 2 - b 2 + b\ <^ n < -b. (21b) 

Since regions ofS! > —6, Vl < —b both exist on 
f2, —b, ft > 0, this direction does not yield any additional 
restrictions. In the reverse direction, we can drop 
the absolute value signs from (21a) since the larger 
quantity remains positive when h > 0. (21b) also 
implies (19) if f2 2 — b 2 + b < 0. Since f2 < — b implies 
il 2 — b 2 + b < for b < 0, (18c) exists everywhere its dis- 
criminant is nonnegative (on the quadrant f2, —b > 0). □ 

Claim: On Q > 0, (18d) exists if and only if ft > 
and SI < —b. 



C. Stability of the fixed points 

At (18a) and (18b), the linearization of (17) has the 
determinant and trace: 

A = b(b+ 1)(1 - n 2 /b 2 ), 

, — — (24) 

T± = ±(26+ \)yj\ -n 2 /b 2 . 

where (18a) has trace r+ and (18b) has trace r_. By 
careful consideration of (24) and Fig. 2, we find that A 
and t+ take signs according to the four cases in Tabic II. 
The case of A and r_ is analogous. 

Similarly, the determinant and trace of the lineariza- 
tion at (18c) and (18d) are 

A ± =±-^-Vh--^-h, (25a) 
1—6 1—6 

r = 0. (25b) 
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where (18c) has determinant A + and (18d) has determi- 
nant A_. We now consider how A± in (25a) takes signs 
as a function of b and SI. The results are summarized in 
Table III. 

Claim: On 17 > 0, (18c) is a center. 
Proof: Clearly, A + > for b < 0. On b > 1, 

161 SI r- 

A+ > ° ^ |1^5| ^ > jl - 6| ■ ^ 

Likewise, for & on (0, 1), 

A+ > °^ ]T^b\^ > \i^~b\ h ^ n > L (27) 

Since SI > b for all (b, SI) where (18c) exists, A + > and 
(18c) is a center on b > 0, as well. □ 

Claim: On SI > 0, (18d) is a center for b > — 1 and a 
saddle for b < — 1. 

Proof: We need only be concerned with the negative 
6-axis, since (18d) docs not exist where SI, fo > 0. For b 
on (-1,0), 

^->0^-J^h>-^—Vh^n<-b, (28) 
|1-&| |1-6| 

while for b < — 1, 

A_ < o- t-^ttV^ > t- ^77 h & SI < -b. (29) 
|1-6| |1-6| 

Since fi < -b for all (b, SI) where (18d) exists, A_ > 
and (18d) is a center for b on (—1,0), while A_ < and 
(18d) is a saddle for b < -1. □ 



D. Parameter space and phase portraits 

Figure 3 summarizes our findings regarding the regions 
of the &-S7 parameter plane where each fixed point exists, 
as well as the stability classifications of the fixed points 
on these regions. If we let SI — > —SI, then |x*| of (18c) 
becomes \x*\ of (18d) and vice versa. Hence, the exis- 
tence and stability of (18c) at a given point (6, —SI) of 
parameter space for SI > is given by the existence and 
stability of (18d) at (b, SI), and vice versa. 

Figure 4 combines the four separate panels of Fig. 3 
into a single image. The various regions in Fig. 3 yield 
six qualitatively distinct regions of the 6- SI plane in Fig. 4. 
We distinguish between the regions (b) and (e) in Fig. 4, 
because the x* at which the single center of these regions 
is located changes sign upon crossing b = — 1. 

Figure 5 plots the phase portraits for the order param- 
eter dynamics governed by (16). The panels show the 
qualitatively different behavior that occurs in the six re- 
gions of Fig. 4. Figure 5(a) depicts what happens in the 



region where b > and SI < b. There, all trajectories 
are attracted to a stable fixed point on the unit circle, 
representing a synchronized rest state. Notice that an in- 
variant vertical line seems to join the repelling fixed point 
with the attracting one. To prove that this vertical line 
truly is invariant, observe from (17) that x — whenever 
x = fl/b. Hence a solution that starts on this line stays 
there forever. 

Figure 5(b) shows the case where b > — 1 and SI > 
The fixed points that previously existed on the unit cir- 
cle have disappeared. They annihilated each other when 
SI = |6|, thereby creating the periodic orbit on the unit 
circle seen in Fig. 5(b). Physically, this orbit represents a 
synchronized oscillation with all the junctions moving in 
phase. But this synchronous state is not attracting; it is 
neutrally stable. In fact, the entire unit disk is filled with 
neutrally stable periodic orbits, all of which surround a 
neutrally stable fixed point, the splay state mentioned 
earlier. 

The remaining panels show examples of additional 
cases when b < 0. We do not dwell on these, as they cor- 
respond to a negative resistance in either the load or the 
Josephson junctions and hence arc physically unrealistic. 
The main features to observe are the saddle connection 
and the coexistence of two neutrally stable splay states 
in Fig. 5(c), and the saddle and center splay states in 
Fig. 5(f). 



IV. EiC-LOADED CIRCUIT 

We turn now to a more complicated kind of Joseph- 
son junction array. Instead of the purely resistive load 
assumed earlier, we allow a load comprised of a resistor, 
inductor, and capacitor in series. This load is placed in 
parallel with TV identical, overdamped Josephson junc- 
tions wired in series. The whole circuit is driven by a 
constant current source lb, as shown in Fig. 6. 

Consider the infinite- TV limit of this system. As shown 
by Strogatz and Mirollo [39] , the time-dependent dynam- 
ics of the array can be written in dimcnsionlcss form as 

v^,t)=- l -e- i * + {I h -Q)+ l -e i * (30) 

where the state variable Q is governed by a dimcnsionlcss 
version of Kirchhoff 's voltage law: 

LQ + {R + l)Q + C- 1 Q = I b - / p(0,t)sin<W (31) 

Jo 

In (30) and (31), Q(t) is the dimensionless charge on the 
capacitor, while the other new quantities are as indicated 
in Fig. 6. 

Equation (30) has the special trigonometric form 
v(4>, t) = je~ 1 ^ + g + fe 1 ^ required by the Ott-Antonsen 
method and therefore is reducible by the method of 
Sec. II. By reading off the / and g implied by (30) and 
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saddle 




FIG. 3: Regions of existence in the upper-half b-£l plane for each of the four fixed points. (A) corresponds to fixed point 
(18a), (B) to fixed point (18b), (C) to (18c) and (D) to (18d). The existence and stability of (18a) at a given point (b, —Q) of 
parameter space for SI > is given by the existence and stability of (18a) at (b, Q), and likewise for (18b). By contrast, the 
existence and stability of (18c) at (6, — Q) is equivalent to the existence and stability of (18d) at (6, fi), and vice versa. 



Q 




FIG. 4: The six qualitatively distinct regions of the b-S} pa- 
rameter plane. The partition is symmetric about the fe-axis. 



substituting them into the Riccati equation (12) for the 
order parameter z, we obtain 



which has real and imaginary parts 



1 



r + r 



■ cos-0 



(33) 



■ sin ip + I b 



By defining P = Q and computing the integral, (31) can 
also be split into a two-dimensional system: 



LP = -(R+1)P - C-h 
Q = P. 



+ lb — r sin ip 



(34) 



Now let P' = LP, Q' = LQ, S = -(R + 1)/L, and 
-0' = 7r/2 — ijj. If we substitute these definitions into (33) 
and (34) and drop the primes, we obtain 



1-r 1 



■ sin?/; 
-1 

-cosip - I b + L~ X P 



P = I b + SP - (LC)- X Q -rcosip 
Q = P. 



(35) 



+ i{h - Q)z, 



This is the low-dimensional system that governs the 

(32) 

flow on the Poisson submanifold. Observe that the state 
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FIG. 5: Representative phase portraits for the six regions of the upper-half b-Sl plane with qualitatively distinct phase plane 
behavior. The trajectories are plotted in polar coordinates r and ip on the unit disk. Solid dots denote Lyapunov stable fixed 
points, while open dots denote unstable fixed points. The letter labels (a)-(f) match with the labels in Fig. 4. 



variables of (35) are the order parameter amplitude r and 
phase ip, along with the dimcnsionlcss rescaled current P 
and charge Q through the load capacitor. The control 
parameters are the bias current Ib, and the load induc- 
tance L, (rescaled) resistance S, and capacitance C . 

A complete analysis of (35) is beyond the scope of 
this paper. From previous numerical experiments, we 
know that there would be many attractors and other 
complicated features to consider [39, 40, 41, 42]. Rather 
than try to enumerate and analyze all of these, our aim 
will be to merely list a few basic facts about the fixed 
points of the system. In particular, we show that an ear- 
lier result — an eigenvalue equation whose roots give the 
four non-trivial Floquet multipliers of the splay state — 
has a more straightforward derivation within the present 
framework. 



A. Fixed points 

We first convert r and ij) of (35) to Cartesian coordi- 
nates x and y as in the purely resistive case. The result 



is 

x = -xy+(I b -L- 1 P)y 

y = \{x 2 -y 2 + l)-{I b -L- l P)x 

2 (36) 
P = I b + SP- {LC)~ l Q-x 

Q = P. 

There are four fixed points of (36): 

x* = I b ; y* = w_; P* = 0; {LC^Q* = (37a) 

x* = I b ; y* = P* = 0; (LC^Q* = (37b) 

x* =I b + u+; y*,P* = 0; Q* = lu+LC (37c) 

x* = I b - uj+; y*,P* = 0; Q* = -uj+LC (37d) 

where uj± = ±I 2 =p 1. Once again, (37a) and (37b) rep- 
resent synchronous fixed points in which the Joscphson 
junctions are in phase and not oscillating, while (37c) and 
(37d) represent splay-state periodic orbits. The deriva- 
tion of (37) makes no assumptions about the parameters 
except that they are real scalars, so we have no need to 
disregard certain parameter values as we did for b in (18). 
Note that (37a) and (37b) exist if and only \I b \ < 1. 
Similarly, the requirement \x*\ < 1 implies (37c) only 
exists on (-co, —1] and (37d) only exists on [l,oo). 
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FIG. 6: A series of Josephson junctions in parallel with an 
RLC load, h represents a dimensionless version of the con- 
stant bias current, C a dimensionless load capacitance, L a 
dimensionless load inductance, and R a dimensionless load 
resistance. 



B. Stability of the splay state 

Since the splay states are the primary concern in much 
of the existing literature, we compute the linearization of 
(36) at (37c): 



( Tw + \ 

±uj + L~ 1 (I b ±u) + ) 

-1 s -(LC)- 1 

\ o o i o / 



J = 



(38) 



The Jacobian (38) has the characteristic polynomial: 

LX 4 + {R+ 1)A 3 + {C- 1 + Llu 2 + )X 2 

+ (I b uj + + Rlu 2 + )X + lj 2 + C~ 1 = (39) 

where we have substituted back in the definition of S and 
multiplied through by L. 

The characteristic polynomial (39) was first derived 
fifteen years ago (see (13) of Rcf. [39]). At the time, its 
derivation gave the first explanation for why there are 
just four non-neutral Floquet multipliers for the splay 
state of the RLC system [41]. It also allowed analytical 
predictions of those multipliers [39]. However, the earlier 
derivation [39] involved Fourier expansions of infinites- 
imal perturbations about the splay states, a procedure 
more complicated than the one given here. 

In retrospect, we can see now that perturbations tan- 
gent to the Poisson submanifold are precisely those re- 
sponsible for the non-neutral directions; perturbations 
transverse to this manifold arc the neutral ones. 



V. DISCUSSION 

It is important to understand both the successes and 
the limitations of our analysis. The systems we have 
studied comprise a special class of Josephson junction 
arrays, namely those in which all the junctions are iden- 
tical and heavily overdamped, meaning that we can ig- 
nore their internal capacitance. The junctions are con- 
nected in series, driven by a constant bias current and 
coupled through a load in parallel. For this class of ar- 
rays, it has been known since 1994 that the governing 
equations specify a family of low-dimensional invariant 
manifolds [44]. For resistively loaded arrays, these in- 
variant manifolds are three-dimensional, while for arrays 
with an RLC load, they are five-dimensional. These re- 
sults hold for any number of junctions and extend to 
infinite N. 

In this paper, we have shown that in the infinite-iV 
limit, the equations can be reduced even more dramati- 
cally. In other words, a strictly smaller invariant subman- 
ifold exists as a degenerate case of the invariant manifolds 
known previously. We have called it the Poisson subman- 
ifold; it consists of all phase distributions taking the form 
of a Poisson kernel. Using the Ott-Antonsen ansatz, we 
explicitly calculated the flow equations on this manifold, 
demonstrating in the process that the Poisson submani- 
fold is invariant. 

The resulting low-dimensional dynamical systems shed 
new light on the behavior of Josephson junction series ar- 
rays. For example, earlier local arguments [38, 39, 40, 41, 
42] showed that the synchronous periodic state and splay 
states for arrays with a resistive load exhibit neutral sta- 
bility to linear order over a wide range of parameters. 
But do they exhibit neutral stability in nonlinear real- 
ity? This has been a long-standing open question. We 
can now observe that the answer is yes. Figure 5(b) shows 
that the splay state is surrounded by neutrally stable pe- 
riodic orbits that fill the Poisson submanifold, continuing 
all the way out to the synchronized orbit on the boundary 
of the disk. 

The other phase portraits in Fig. 5 similarly provide 
new information regarding the global structure. For ex- 
ample, they elucidate how splay states bifurcate with one 
another and with the synchronous periodic state and re- 
veal global features such as the vertical hcteroclinic or- 
bit joining the in-phase rest states in several panels of 
Figure 5. Previous global results regarding Josephson 
junction arrays were confined to the averaged versions of 
such systems [43, 44], which were derived via perturba- 
tion methods in the limit of weak coupling or high bias 
current [49]. Our results, by contrast, hold for all val- 
ues of the circuit parameters. The trade-off is that they 
require infinite N. 

However, the most serious drawback of our analysis 
is that it focuses on a thin slice of phase space which 
is unrepresentative of the dynamics of the full system. 
For instance, the Poisson submanifold for the resistively 
loaded array is a degenerate, two-dimensional leaf in the 
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foliation of phase space by three-dimensional invariant 
manifolds. Could attractors or other interesting dynam- 
ical states exist off the Poisson submanifold? Numerical 
simulations say yes: KAM-like chaos occurs in the origi- 
nal equations for resistively loaded arrays [42, 44]. 

The remaining challenge is then to show that the 
reduced equations faithfully capture this chaos on the 
larger invariant manifolds. We will address this issue in 



a subsequent paper [48] in which we use Mobius transfor- 
mations to simplify the circuit equations for Joscphson 
junction arrays. Essentially the same idea has been de- 
veloped independently by Pikovsky and Rosenblum [47], 
who have obtained new results for the Kuramoto model 
as well as more complex hierarchies of oscillators. 

Acknowledgments: Research supported in part by 
National Science Foundation grant NSF CISE-0835706. 



A. T. Winfree, J. Theor. Biol. 16, 15 (1967). [26] 
Y. Kuramoto, Chemical Oscillations, Waves, and Turbu- 
lence (Springer, Berlin, 1984). [27 
J. A. Acebron, L. L. Bonilla, C. J. P. Vicente, F. Ritort, [28 
and R. Spigler, J. Theor. Biol. 77, 137 (2005). 
S. H. Strogatz, in Frontiers in Mathematical Biology, [29 
edited by S. A. Levin (Springer, New York, 1994), vol. 
100 of Lecture Notes in Biomathematics, pp. 122-138. 
N. Wiener, Nonlinear Problems in Random Theory (MIT [30 
Press, Cambridge, MA, 1958). 

N. Wiener, Cybernetics (MIT Press, Cambridge, MA, [31 
1961), 2nd ed. 

T. D. Frank, A. Daffertshofera, C. E. Pepera, P. J. Beeka, [32 
and H. Hakenb, Physica D 144, 62 (2000). 

C. S. Peskin, Mathematical Aspects of Heart Physiology [33 
(Courant Institute of Mathematical Sciences, New York, 
1975). [34 

D. C. Michaels, E. P. Matyas, and J. Jalife, Circ. Res. [35 
58, 706 (1986). 

D. C. Michaels, E. P. Matyas, and J. Jalife, Circ. Res. [36 
61, 704 (1987). 

J. Buck and E. Buck, Science 159, 1319 (1968). [37 

F. E. Hanson, Fed. Proc. 37, 2158 (1978). 
J. Buck, Quart. Rev. Biol. 63, 265 (1988). [38 

G. B. Ermentrout, J. Math. Biol. 29, 571 (1991). 
D. Kim, BioSystems 76, 7 (2004). [39 
M. K. McClintock, Nature 229, 244 (1971). 
M. J. Russell, G. M. Switz, and K. Thompson, Pharma- [40 
col. Biochem. Behav. 13, 737 (1980). 

H. C. Wilson, Psychoneuroendocrinology 17, 565 (1992). [41 
A. K. Ghosh, B. Chance, and E. K. Pye, Arch. Biochem. 
Biophys. 145, 319 (1971). [42 
D. Njus, V. D. Gooch, and J. W. Hastings, Cell Biophys. 
3, 223 (1981). [43 
J. W. Hastings, H. Broda, and C. H. Johnson, in Tem- 
poral Order, edited by L. Rensing and N. I. Jaeger [44 
(Springer, Berlin, 1985), pp. 213-221. 

S. Dano, F. Hynne, S. D. Monte, F. d'Ovidio, P. G. [45 
Sorensen, and H. Westerhoff, Faraday Discuss. 120, 261 [46 
(2001). [47 
D. B. Murray, S. Roller, H. Kuriyama, and D. Lloyd, J. 
Bacteriol. 183, 7253 (2001). [48 
J. Garcia-Ojalvo, M. B. Elowitz, and S. H. Strogatz, Pro- 
ceed. Nat. Acad. Sci. 101, 10955 (2004). [49 
H. G. Winful and S. S. Wang, Appl. Phys. Lett. 53, 1894 
(1988). 



K. Wiesenfeld, C. Bracikowski, G. James, and R. Roy, 
Phys. Rev. Lett. 65, 1749 (1990). 

R.-D. Li and T. Erneux, Phys. Rev. A 46, 4252 (1992). 
L. Fabiny, P. Colet, R. Roy, and D. Lenstra, Phys. Rev. 
A 47, 4287 (1993). 

S. Y. Kourtchatov, V. V. Likhanskii, A. P. . Napartovich, 

F. T. Arecchi, and A. Lapucci, Phys. Rev. A 52, 4089 
(1995). 

G. Kozyreff, A. G. Vladimirov, and P. Mandel, Phys. 
Rev. Lett. 85, 3809 (2000). 

T. Heil, I. Fischer, W. Elsasser, J. Mulet, and C. R. Mi- 
rasso, Phys. Rev. Lett. 86, 795 (2001). 
S. H. Strogatz, C. M. Marcus, and R. M. Westervelt, 
Phys. Rev. Lett. 61, 2380 (1988). 

S. H. Strogatz and R. M. Westervelt, Phys. Rev. B 40, 
10501 (1989). 

A. A. Middleton, Phys. Rev. Lett. 68, 670 (1992). 

W. Wang, I. Z. Kiss, and J. L. Hudson, Chaos 10, 248 

(2000). 

I. Z. Kiss, Y. Zhai, and J. L. Hudson, Phys. Rev. Lett. 
88, 238301 (2002). 

I. Z. Kiss, Y. Zhai, and J. L. Hudson, Science 296, 1676 
(2002). 

K. Y. Fsang, R. E. Mirollo, S. H. Strogatz, and 

K. Wiesenfeld, Physica D 48, 102 (1991). 

S. H. Strogatz and R. E. Mirollo, Phys. Rev. E 47, 220 

(1993) . 

K. Y. Tsang and I. B. Schwartz, Phys. Rev. Lett. 68, 
2265 (1992). 

S. Nichols and K. Wiesenfeld, Phys. Rev. A 45, 8430 
(1992). 

D. Golomb, D. Hansel, B. Shraiman, and H. Somopolin- 
sky, Phys. Rev. A 45, 3516 (1992). 

S. Watanabe and S. H. Strogatz, Phys. Rev. Lett. 70, 
2391 (1993). 

S. Watanabe and S. H. Strogatz, Physica D 74, 194 

(1994) . 

C. J. Goebel, Physica D 80, 18 (1995). 

E. Ott and T. M. Antonsen, Chaos 18, 037113 (2008). 
A. Pikovsky and M. Rosenblum, arXiv:0809.3700v2 
[nlin.AO] (2008). 

R. Mirollo, S. A. Marvel, and S. H. Strogatz (2009), forth- 
coming. 

J. W. Swift, S. H. Strogatz, and K. Wiesenfeld, Physica 
D 55, 239 (1992). 



